home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / rayshade / libray / libobj / roots.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  5KB  |  249 lines
  1. /*
  2.  *  Roots3And4.c
  3.  *
  4.  *  Utility functions to find cubic and quartic roots,
  5.  *  coefficients are passed like this:
  6.  *
  7.  *      c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0
  8.  *
  9.  *  The functions return the number of non-complex roots and
  10.  *  put the values into the s array.
  11.  *
  12.  *  Author:         Jochen Schwarze (schwarze@isa.de)
  13.  *
  14.  *  Jan 26, 1990    Version for Graphics Gems
  15.  *  Oct 11, 1990    Fixed sign problem for negative q's in SolveQuartic
  16.  *                  (reported by Mark Podlipec),
  17.  *                  Old-style function definitions,
  18.  *                  IsZero() as a macro
  19.  *  Nov 23, 1990    Some systems do not declare acos() and cbrt() in
  20.  *                  <math.h>, though the functions exist in the library.
  21.  *                  If large coefficients are used, EQN_EPS should be
  22.  *                  reduced considerably (e.g. to 1E-30), results will be
  23.  *                  correct but multiple roots might be reported more
  24.  *                  than once.
  25.  */
  26.  
  27. #include "libcommon/common.h"
  28.  
  29. extern double   sqrt(), cbrt(), cos(), acos();
  30.  
  31. /* epsilon surrounding for near zero values */
  32.  
  33. /*
  34.  * In case M_PI isn't defined in math.h
  35.  */
  36. #ifndef M_PI
  37. #define M_PI PI
  38. #endif
  39.  
  40. #define     EQN_EPS     1e-9
  41. #define        IsZero(x)    ((x) > -EQN_EPS && (x) < EQN_EPS)
  42.  
  43. #ifndef CBRT
  44. #define     cbrt(x)     ((x) > 0.0 ? pow((double)(x), 1.0/3.0) : \
  45.               ((x) < 0.0 ? -pow((double)-(x), 1.0/3.0) : 0.0))
  46. #endif
  47.  
  48.  
  49. int SolveQuadric(c, s)
  50.     double c[ 3 ];
  51.     double s[ 2 ];
  52. {
  53.     double p, q, D;
  54.  
  55.     /* normal form: x^2 + px + q = 0 */
  56.  
  57.     p = c[ 1 ] / (2 * c[ 2 ]);
  58.     q = c[ 0 ] / c[ 2 ];
  59.  
  60.     D = p * p - q;
  61.  
  62.     if (IsZero(D))
  63.     {
  64.     s[ 0 ] = - p;
  65.     return 1;
  66.     }
  67.     else if (D > 0)
  68.     {
  69.     double sqrt_D = sqrt(D);
  70.  
  71.     s[ 0 ] =   sqrt_D - p;
  72.     s[ 1 ] = - sqrt_D - p;
  73.     return 2;
  74.     }
  75.     else /* if (D < 0) */
  76.         return 0;
  77. }
  78.  
  79.  
  80. int SolveCubic(c, s)
  81.     double c[ 4 ];
  82.     double s[ 3 ];
  83. {
  84.     int     i, num;
  85.     double  sub;
  86.     double  A, B, C;
  87.     double  sq_A, p, q;
  88.     double  cb_p, D;
  89.  
  90.     /* normal form: x^3 + Ax^2 + Bx + C = 0 */
  91.  
  92.     A = c[ 2 ] / c[ 3 ];
  93.     B = c[ 1 ] / c[ 3 ];
  94.     C = c[ 0 ] / c[ 3 ];
  95.  
  96.     /*  substitute x = y - A/3 to eliminate quadric term:
  97.     x^3 +px + q = 0 */
  98.  
  99.     sq_A = A * A;
  100.     p = 1.0/3 * (- 1.0/3 * sq_A + B);
  101.     q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
  102.  
  103.     /* use Cardano's formula */
  104.  
  105.     cb_p = p * p * p;
  106.     D = q * q + cb_p;
  107.  
  108.     if (IsZero(D))
  109.     {
  110.     if (IsZero(q)) /* one triple solution */
  111.     {
  112.         s[ 0 ] = 0;
  113.         num = 1;
  114.     }
  115.     else /* one single and one double solution */
  116.     {
  117.         double u = cbrt(-q);
  118.         s[ 0 ] = 2 * u;
  119.         s[ 1 ] = - u;
  120.         num = 2;
  121.     }
  122.     }
  123.     else if (D < 0) /* Casus irreducibilis: three real solutions */
  124.     {
  125.     double phi = 1.0/3 * acos(-q / sqrt(-cb_p));
  126.     double t = 2 * sqrt(-p);
  127.  
  128.     s[ 0 ] =   t * cos(phi);
  129.     s[ 1 ] = - t * cos(phi + M_PI / 3);
  130.     s[ 2 ] = - t * cos(phi - M_PI / 3);
  131.     num = 3;
  132.     }
  133.     else /* one real solution */
  134.     {
  135.     double sqrt_D = sqrt(D);
  136.     double u = cbrt(sqrt_D - q);
  137.     double v = - cbrt(sqrt_D + q);
  138.  
  139.     s[ 0 ] = u + v;
  140.     num = 1;
  141.     }
  142.  
  143.     /* resubstitute */
  144.  
  145.     sub = 1.0/3 * A;
  146.  
  147.     for (i = 0; i < num; ++i)
  148.     s[ i ] -= sub;
  149.  
  150.     return num;
  151. }
  152.  
  153.  
  154. int SolveQuartic(c, s)
  155.     double c[ 5 ]; 
  156.     double s[ 4 ];
  157. {
  158.     double  coeffs[ 4 ];
  159.     double  z, u, v, sub;
  160.     double  A, B, C, D;
  161.     double  sq_A, p, q, r;
  162.     int     i, num;
  163.  
  164.     /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
  165.  
  166.     A = c[ 3 ] / c[ 4 ];
  167.     B = c[ 2 ] / c[ 4 ];
  168.     C = c[ 1 ] / c[ 4 ];
  169.     D = c[ 0 ] / c[ 4 ];
  170.  
  171.     /*  substitute x = y - A/4 to eliminate cubic term:
  172.     x^4 + px^2 + qx + r = 0 */
  173.  
  174.     sq_A = A * A;
  175.     p = - 3.0/8 * sq_A + B;
  176.     q = 1.0/8 * sq_A * A - 1.0/2 * A * B + C;
  177.     r = - 3.0/256*sq_A*sq_A + 1.0/16*sq_A*B - 1.0/4*A*C + D;
  178.  
  179.     if (IsZero(r))
  180.     {
  181.     /* no absolute term: y(y^3 + py + q) = 0 */
  182.  
  183.     coeffs[ 0 ] = q;
  184.     coeffs[ 1 ] = p;
  185.     coeffs[ 2 ] = 0;
  186.     coeffs[ 3 ] = 1;
  187.  
  188.     num = SolveCubic(coeffs, s);
  189.  
  190.     s[ num++ ] = 0;
  191.     }
  192.     else
  193.     {
  194.     /* solve the resolvent cubic ... */
  195.  
  196.     coeffs[ 0 ] = 1.0/2 * r * p - 1.0/8 * q * q;
  197.     coeffs[ 1 ] = - r;
  198.     coeffs[ 2 ] = - 1.0/2 * p;
  199.     coeffs[ 3 ] = 1;
  200.  
  201.     (void) SolveCubic(coeffs, s);
  202.  
  203.     /* ... and take the one real solution ... */
  204.  
  205.     z = s[ 0 ];
  206.  
  207.     /* ... to build two quadric equations */
  208.  
  209.     u = z * z - r;
  210.     v = 2 * z - p;
  211.  
  212.     if (IsZero(u))
  213.         u = 0;
  214.     else if (u > 0)
  215.         u = sqrt(u);
  216.     else
  217.         return 0;
  218.  
  219.     if (IsZero(v))
  220.         v = 0;
  221.     else if (v > 0)
  222.         v = sqrt(v);
  223.     else
  224.         return 0;
  225.  
  226.     coeffs[ 0 ] = z - u;
  227.     coeffs[ 1 ] = q < 0 ? -v : v;
  228.     coeffs[ 2 ] = 1;
  229.  
  230.     num = SolveQuadric(coeffs, s);
  231.  
  232.     coeffs[ 0 ]= z + u;
  233.     coeffs[ 1 ] = q < 0 ? v : -v;
  234.     coeffs[ 2 ] = 1;
  235.  
  236.     num += SolveQuadric(coeffs, s + num);
  237.     }
  238.  
  239.     /* resubstitute */
  240.  
  241.     sub = 1.0/4 * A;
  242.  
  243.     for (i = 0; i < num; ++i)
  244.     s[ i ] -= sub;
  245.  
  246.     return num;
  247. }
  248.  
  249.